/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"

#include "uhm/matrix/uhm/matrix.hxx"
#include "uhm/matrix/uhm/fla.hxx"

namespace uhm {
  void Matrix_FLA_::check_ldl_nopiv_1() {
#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix 
    // ----------------------------------------------------------
    if (this->fs) {
      FLA_Copy( ~(this->flat.xt), ~(this->flat.rt) );

      FLA_Hermitianize( FLA_LOWER_TRIANGULAR, ~(this->flat.ATL) );
      linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, 
                                              FLA_UPPER_TRIANGULAR, ~(this->flat.ATL) );
      linal::trmm( FLA_LEFT, FLA_UPPER_TRIANGULAR, FLA_NO_TRANSPOSE,
                          FLA_NONUNIT_DIAG, FLA_ONE, 
                          this->hier.ATL, this->hier.rt ); 
    }

    if (this->fs && this->ss) {
      int bmn = get_hier_blocksize();
      
      this->flat.ATR.create( this->flat.ABL.get_datatype(), 
                             this->flat.ABL.get_n(),
                             this->flat.ABL.get_m() );
      this->hier.ATR.create( this->flat.ATR, bmn, bmn );
      FLA_Copyt( FLA_TRANSPOSE, ~(this->flat.ABL), ~(this->flat.ATR) );
      linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, 
                                              ~(this->flat.ATL), ~(this->flat.ATR) );
      
      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, FLA_ONE,
                          this->hier.ATR, this->hier.xb, 
                          FLA_ONE, this->hier.rt );

      this->hier.ATR.free();
      this->flat.ATR.free();
    }
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    if (this->fs) {
      FLA_Copy( ~(this->flat.xt), ~(this->flat.rt) );

      FLA_Hermitianize( FLA_LOWER_TRIANGULAR, ~(this->flat.ATL) );
      linal::apply_diag_scal( FLA_LEFT, FLA_LOWER_TRIANGULAR, this->flat.ATL );
      FLA_Trmm( FLA_LEFT, FLA_UPPER_TRIANGULAR, FLA_NO_TRANSPOSE,
                FLA_NONUNIT_DIAG, FLA_ONE, 
                ~(this->flat.ATL), ~(this->flat.rt) ); 
    }

    if (this->fs && this->ss) {
      this->flat.ATR.create( this->flat.ABL.get_datatype(), 
                             this->flat.ABL.get_n(),
                             this->flat.ABL.get_m() );
      FLA_Copyt( FLA_TRANSPOSE, ~(this->flat.ABL), ~(this->flat.ATR) );
      linal::apply_diag_scal( FLA_LEFT, this->flat.ATL, this->flat.ATR );
      
      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, FLA_ONE,
                ~(this->flat.ATR), ~(this->flat.xb), 
                FLA_ONE, ~(this->flat.rt) );

      this->flat.ATR.free();
    }

#endif
  }
   
  // from leaf to root
  void Matrix_FLA_::check_ldl_nopiv_2() {
#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix 
    // ----------------------------------------------------------
    if (this->fs && this->ss) {
      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, FLA_ONE,
                          this->hier.ABL, this->hier.rt, 
                          FLA_ONE, this->hier.rb );
    }
    
    if (this->fs) {
      linal::trmm( FLA_LEFT, FLA_LOWER_TRIANGULAR, FLA_NO_TRANSPOSE,
                          FLA_UNIT_DIAG, FLA_ONE, 
                          this->hier.ATL, this->hier.rt ); 
    }
    // rb should be merged for upper hierarchy
    // pivot should be applied before it is merged
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    if (this->fs && this->ss) {
      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, FLA_ONE,
		~(this->flat.ABL), ~(this->flat.rt), 
		FLA_ONE, ~(this->flat.rb) );
    }
    
    if (this->fs) {
      FLA_Trmm( FLA_LEFT, FLA_LOWER_TRIANGULAR, FLA_NO_TRANSPOSE,
		FLA_UNIT_DIAG, FLA_ONE, 
		~(this->flat.ATL), ~(this->flat.rt) ); 
    }
    // rb should be merged for upper hierarchy
    // pivot should be applied before it is merged
#endif
  }
}
